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Abstract 

The limit cycle of the van der Pol oscillator, x + e(x 2 — l)x + x = 0, is studied 
in the plane (x, x) by applying the homotopy analysis method. A recursive set of 
formulas that approximate the amplitude and form of this limit cycle for the whole 
range of the parameter e is obtained. These formulas generate the amplitude with 
an error less than 0.1%. To our knowledge, this is the first time where an analytical 
approximation of the amplitude of the van der Pol limit cycle, with validity from 
the weakly up to the strongly nonlinear regime, is given. 
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1 Introduction 

The differential equation, 

x{t) + e{x 2 - l)x(t) +x(t) = 0, t>0, (1.1) 

with e a real parameter and the dot denoting the derivative with respect to time 
t, is called the van der Pol oscillator For e > 0, and due to the nonlinear term 
e(x 2 — l)x, the system accumulates energy in the region | x |< 1 and dissipates 
this energy in the region | x |> 1. This constraint implies the existence of an stable 
periodic motion {limit cycle [2j) when e > 0. If the nonlinearity is increased, the 



dynamics in the time domain runs from near-harmonic oscillations when e — > to 
relaxation oscillations when e — > oo, making it a good model for many practical 
situations [3114] . The closed curve representing this oscillation in the plane (x,x) is 
quasi circular when e — ► and a sharp figure when e — > oo. For e < 0, the dynamics 
is dissipative in the region | x |< 1 and amplificative for | x |> 1. Under these 
conditions, the periodic motion is still possible but unstable. In this case, the limit 
cycle can be derived from that one with e > taking into account the symmetry 
(e,x(t)) — ► (— e, —x(—t)). Therefore, it is enough to study the case e > to obtain 
also the behavior of the system for e < 0. 



Different standard methods (perturbative, non-perturbative, geometrical) [1,2,3,4, 5. 6, 
have permitted to study extensively the limit cycle of the van der Pol equation, in 
the weakly (e — > 0) and in the strongly (e — > oo) nonlinear regimes. However, in- 
vestigations giving analytical information of this object in the intermediate regime 
(0 < e C oo) are lacking in the literature. In this paper, it is our aim to fill in 
this gap by applying to the Eq. (1.1) the homotopy analysis method (HAM) intro- 
duced by Liao [8, 9) in the nineties. This method has been shown to be powerful to 
solve different nonlinear problems |10llll|ll2|ll3|ll4|ll5|ll6j . In particular, it has been 
applied in [T7| to Lienard equation, x + ef(x)x + x = 0, which is the generalization 
of the van der Pol system when f(x) is an arbitrary function. As the interest in 
that work [17] was the amplitude and the frequency of the periodic motions, the 
calculations were performed with the time variable being explicit. As here we are 
only interested in the amplitude and form of the limit cycles, the time dependence 
of the solutions can be omitted, and we can work directly in the phase space (x, x). 
Our results for the van der Pol limit cycle are presented in Section 2. In Section 3 
it is explained how these results have been obtained from our specific application 
of the HAM to the van der Pol equation. Some conclusions are given in the last 
Section. 



2 The amplitude of the van der Pol limit cycle 



Nowadays, the amplitude of the van der Pol limit cycle a can be easily approximated 
with a classical Runge-Kutta method. The results of this computational calculation 
cle are shown in Table 1. Let us observe that ag(e) > 2 for e > 0, with the asymptotic 
values: a^(e — > 0) = a^(e — > oo) = 2. An upper bound rigorously established in [18] 
for the amplitude a is 2.3233. However, as it was signalled there [18], the maximum 
of cle is 2.02342 and it is obtained for e = 3.3. In view of this result, Odani [18J 
conjectured that the amplitude of the limit cycle of the van der Pol equation is 
estimated by 2 < a(e) < 2.0235 for every e > 0. 
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Table 1. The value a E represents the amplitude a of the van der Pol limit cycle 



obtained by integrating directly Eq. (1.1) with a Runge-Kutta method for the 
indicated values of e. 

Due to the difficulty of calculating the amplitude a by analytical means, one could 
propose as a good solution the constant amplitude 

a = 2 (2.1) 

for the whole range of the parameter e. Knowing that the experimental upper bound 
for a is 2.02342, i.e. 2 < a(e) < 2.02342 for every e > 0, the error made with this 
approximation, (a(e) — a)/a(e), would be about 1%. 

If more precision is needed, the different analytical expansions in e that have been 
found for the amplitude in the weakly and in the strongly nonlinear regimes can 
be considered. Evidently, the error becomes very large in the regions where these 
approximations are not valid. In [1749 20], a recursive perturbation approximation 
is used to find the formula for the amplitude when e — > 0. This is, up to order 0(e 8 ), 
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(2.2) 



This expansion agrees for small e with the computational calculation a E presented 
in Table 1. For < e < 2 the error is less than 1%, for e ~ 4 the error is bigger 
than 10%, and for e ~ 6 the formula has lost its validity and the error is bigger 
than 50%. In [21], the asymptotic dependence of the amplitude on e is given for 
sufficiently large e > 0, 



a e 



2 + 0.7793 e~ 4/3 + 0(e 



4/3N 



(2.3) 



where O means (2(e~ 4 / 3 ) < (9(e -4 / 3 ). Compared with a E , this formula generates 
the amplitude with an error bigger than 15% for e < 2. The formula starts to have 
validity for e ~ 10 with an error around 1%, that passes to be less than 0.1% when 



e > 50. In Figs, l.a and l.b, formulas (2.2) and (2.3) are respectively plotted versus 
e in their regions of validity. 



We propose here two formulas, a^i(e) and a R2 (e), that approximate the amplitude 
of the van der Pol limit cycle for every e > 0. The first one with an error less than 
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(a) (b) 

Fig. 1. Comparison of the 'experimental' amplitude ae(e) (dotted curve) with (a) am- 
plitude a(e) e ^o given by formula (2.2) (solid line), and (b) amplitude a(e) e _K>o given by 



formula (2.3) (solid line). 

0.1%, and the second one with an error less than 0.05%. These formulas have been 
obtained by applying the HAM to the van der Pol equation. The details of the 
derivation of these formulas are given in Section 3. 



The first formula is 



1.737 e 2 



am{e) 



(87r + 9e)(4 + e 2 



(2.4) 



that derives from expression (3.30). The error obtained, \a(e) — a R1 (e)\/a(e), with 



this formula is less than 0.1% for every e > 0. Fig. 2.1 shows am in comparison with 
the 'experimental' amplitude ag. The maximum of am is 2.02317 and it is taken 
for e = 3.3. 



e 2 (75.3562 + 43.0023e + 28.1589e 2 + 8.3479e 3 
(87r + 9e) 2 (4 + e 2 ) 2 



The second formula is 

0.74958 e 2 
am[e) ~ 2 + (8 7 r + 9e)(4 + e 2) 

(2.5) 

that derives from expression (3.33). The error obtained, \a(e) — a R2 (e)\/a(e), with 
this formula is less than 0.05% for every e > 0. The maximum of am is 2.02346 
and it is taken for e = 3.29482. The plot of am in comparison with the amplitude 
aE obtained computationally can be seen in Fig. 2b. Let us observe that it can be 
guessed from Fig. 2b that a# 2 (e) > a E (e) for every e > 0. 



The expansion of expression (2.5) for small e gives 



o^h) = 2 + 0.01491 e 2 + 0(e 3 ). 



For large e, we obtain 



2 + 0.18648 e- L + C(e 



(2.6) 



(2.7) 



Let us note the different scaling of this last expression (with behavior e r ) respect 



to expansion (2.3) (with behavior e L33 ). Taking into account that these approxi- 
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Fig. 2. (a) Amplitude ORi(e) given by formula (2.4) (solid curve), and experimental ampli- 
tude oe(e). (b) Amplitude am(^) given by formula (2.5) (solid curve), and experimental 
amplitude ae(e) (dotted curve), and e given in logarithmic scale. 

mations start to be valid when e > 100, it can be easily seen that this difference is 
negligible, in fact less than 0.01%. Nevertheless, we have tried to modify formula 



(2.5) to obtain the correct scaling e of expression (2.3), but the increase of the 



error in other regions of the parameter e does not recommend this possibility. 

In summary, let us remark the exceptional fit of the 'experimental' points by 



the amplitudes am and am generated with formulas (2.4) and (2.5), respectively 



(see Figs. 2a and 2b). Moreover, by inspection of Fig. 2b, let us finish this section 
by posing the following 



Conjecture: the amplitude am( e ) obtained with formula (2.5) is an upper bound 
for the amplitude a(e) of the van der Pol limit cycle, that is 

o-m( e ) > a ( e ) / or every e > 0. (2.8) 



3 The HAM and the van der Pol equation 



The generalization of the van der Pol equation, u + e(u 2 — l)u + u = 0, is called the 
Lienard equation. In coordinates (u, v), it reads 

u(t)+ef(u)ii(t)+u(t) = 0, t>0, (3.1) 

where we consider that f{u) is an even function. The HAM has been applied to this 
equation in [T7] working explicitly in the time domain. We proceed now to apply 
the HAM to this system by following a different strategy, namely, by omitting the 



time variable of Eq. (3.1). 



3.1 HAM applied to Lienard equation 



If we define the variable v = ii, and suppose that v is a function of u, the acceleration 
u can be rewritten as v'v, where the prime denotes the derivative respect to u. The 
result is a new equation with the time variable eliminated: 

vv' + ef(u)v + u = 0. (3.2) 

When f(u) is even [221123] . the limit cycles of this equation are closed trajectories 
of amplitude a in the (u,v) plane, with —a <u<a and v(a) = v(—a) = 0. 

After the re-scaling (u, v) = (ax, ay) to the new variables (x, y), the limit cycles are 
transformed in solutions of the equation 

yy' + eyf(ax)+x = 0, -1 < x < 1, y{-\) = y(l) = 0, (3.3) 

where the prime is now denoting the derivative respect to x. 

From [22.23J, we know that an approximate solution of the positive branch solution 



of (3.3), for any value of e > 0, is 



, e if — 1 < x < xnl a 

y{x) = Vl - x 2 + - { - 1 (3.4) 

a [ F(a) - F(ax) if x /a < x < 1, 

where F\x) = f(x), x is a negative root of f(x): f(x ) = with x < 0, and a is 
a positive root of F(x ) — F(x): F(x ) = F(a) with a > 0. The number of possible 
stable limit cycles in the strongly non-linear regime (for large e) is obtained from 
the number of positive roots of F(xo) — F(x) when xq runs over the set of all the 
negative roots of f(xo) = 0. We take the positive roots a selected by the algorithm 
explained in [2"2"f2"3"] as the corresponding approximate amplitudes. Let us observe 
that the initial approximate limit cycle y(x) recovers the exact circular form when 
e — >• 0, and also, when e — > oo, it becomes exactly the two-piecewise limit cycle 
proposed in [22]. Moreover, it must be noticed that y(x) has sense only when e > 0, 
and then the results obtained from this initial guess only will have validity for e > 0. 

After eliminating the time variable, the non-linear differential equation yy'+eyf(ax)+ 
x = is of order one. Then, to construct the homotopy [8|9] , we build a linear differ- 
ential equation of order one for the whole range [0, 1] of what is called the homotopy 
parameter p. 

We define an auxiliary linear operator £ by 

d 

C[<p(x]p)} = —<p(x]p), (3.5) 

with the property 

C[d]=Q, (3.6) 



where C\ is a constant, and p is a (homotopy) parameter explained below. 



From (3.3) we define a nonlinear operator 

N[<j>{x-p),A{p)] = <p( X] p) d ^^ + e <j ) (x ] p)f(A(p)x) + x, (3.7) 

and then construct the homotopy 

H[<t>{x;p),A(p)\ = (l-p)£[<j>(x;p)-y(x)] + hpN[(j>{x;p), A{p)}, (3.8) 

where h is a nonzero auxiliary parameter. Setting Tt[(f)(x;p), A(p)} = 0, we have the 
zero-order deformation equation 

(1 -p)C[4>{x-p) - y{x)) + hpM[(f>{x;p),A{p)] = 0, (3.9) 

subject to the boundary conditions 

0(-l;p) = O, 0(l;p) = O, (3.10) 

where p E [0, 1] is an embedding parameter. When the parameter p increases from 
to 1, the solution 4>(x;p) varies from yo(x) = y(x) to y(x), A{p) varies from a to 
a. Assume that <p(x;p) and A(p) are analytic in p e [0, 1] and can be expanded in 
the Maclaurin series of p as follows: 

+oo +oo 

(f>{x;p) = £ y n {x)p n , A(p) = ^ a n p n , (3.11) 

n=0 n=0 



where 

. . 1 d n (j)(x;p) 



n\ dp n 



1 d n A(p) 
p=q n n\ dp n 



p=0 



Notice that series (3.11) contain the auxiliary parameter h, which has influence on 
their convergence regions. Assume that h is properly chosen such that all of these 
Maclaurin series are convergent at p = 1. Hence at p = 1 we have 

+oo +oo 

y(x) = y (x) + Vn(x), a = a + J2 a n- 

n=l n=l 



At the A^th-order approximation, we have the analytic solution of Eq. (3.3), namely 

N N 

y(x) w Y N (x) := ^ Vn{x), a^A N :=J2 a n- (3.12) 

n=0 n=0 



The parameter h is free and can be chosen arbitrarily, in particular, it can be a 
function of e. Nevertheless, the function y{x) is the exact limit cycle solution of 
the Lienard equation (3.3) in the limits e — » and e — ► oo. This means that the 
general solution y(x,h) of Eq. (3.9) should tend to y(x) for any p in those limits 
of e, loosing its dependence on h. Then, a reasonable property for h would be to 



vanish in those limits of e. Hence, the solution of (3.9) would be the exact solution 
y(x) for any value of the parameter p in the limits e — * and e 
function h(e) satisfying these properties is 



oo. A simple 



h(e) 



be 
c + e 



2 - 



with b.c G R. 



(3.13) 



Differentiating Eqs. (3.9) and (3.10) n times with respect to p, then setting p — 0, 
and finally dividing by n\ , we obtain the nth-order deformation equation 



0. 



n 



C[y n (x) - XnVn-iix)] + hR n (x 
subject to the boundary conditions 

y n {-l) = 0, y n {l) = 0, 

where R n (x) is defined by 

1 d n -W[(p(x;p),A(p)} 



1,2,3, 



RJx) 



[n 



dp 



,n— 1 



p=0 



(3.14) 



(3.15) 



(3.16) 



and 



Xn 




At each iteration, we have two unknowns, Ci, in Eq. (3.6), and a n . These unknowns 
are obtained by considering the boundary conditions (3.15) as follows. 

At zero order, we obtain 



yo(x) 



y(x). 



At lth-order, we obtain 



yi{x) 



h 



x 



Vo( x ) 



he / y (t)f(a t)dt, 



-i 



(3.17) 



(3.18) 



where we have imposed the condition yi(— 1) = choosing the integration constant, 
Ci, appropriately. At this moment ao is free, but we can fix the value of by 
imposing y\(l) = 0: 

(3.19) 



J yo(t)f(aot)dt = 0. 



This is a non linear equation for ao- When f(x) is an even polynomial of degree 2n, 
it is a polynomial equation of degree n in the variable Oq. Therefore, at this order, 
the number of limit cycles is the number of positive roots of this equation, at most 
n (see [23] for a longer discussion of this point). 



Then, for every one of the ao solutions of the above equation, i.e. for every one 
of the limit cycles of the system, we proceed order by order in p to obtain higher 
order corrections to the amplitude a and to the shape of the limit cycle y{x). At 



any nth-order, with n > 2, we obtain the nth-order correction y n to the shape of 
the limit cycle: 



h n ~ 1 r x 
y n (x) = y n -i(x) - - V y k (x)y n - k -i{x) -eh y n -i(t)f(a t)dt 

2 k% J - 1 



~ eh H / Vn-k-l(t) 

k=l J ~ 1 m=l 



B ktrn (a 1 ,a 2 , ...,a k j m , (m) 



(3.20) 



m! 



t m f (m) (a t)dt, 



where 5„ jm (ai, a n ) are the partial ordinary Bell polynomials [[21], p. 190]. We 
have £? ,o = 1 an d, for n = 1, 2, 3, B n+1)0 = and they satisfy the recurrence 



n-l 



(ai,...,a n )— 2J a n-kBk ATI— 1 (ai, ... 

k=m— 1 



(3.21) 



We recall that -B n , m (ai, ...,a n ) is linear in a n . 

At nth-order, we can obtain also the nth-order correction a n _i to the amplitude of 
the limit cycle by imposing the condition y n (l) = 0. The result is a linear equation 
for a n _i with solution: 



a n -i 



/X n— 1 fc * 

y n -i{x)f(a x)dx + Vn-k-i 

' 1 k=l m=l 



X 



ml 



x m f m) (a x)dx 



-i 



xy Q (x)f'(a x)dx 



(3.22) 

where the star in the sum symbol means that the term (k,m) = (n — 1,1) is 
excluded. 

For example, at 2th-order we obtain 



SfeCaO = 2/1 (a;) - hy {x)yi(x) - he J [yx(t)f(a t) + ty^fXa^a^dt (3.23) 



and 



ai 



yi(x)f(a x)dx 



xy {x)f'(a x)dx 



(3.24) 



Therefore, a first approximation to the shape of the limit cycle is yo(x) +y 1 (x), and 
a first order approximation to the amplitude of the limit cycle is a — ao + o,\- 

At 3th-order, we have 



Vz{x) = y 2 (x) - hy (x)y 2 (x) - ^yl(x) 



-he 



1 



y2(t)f(a Q t) + aityi(t)f'(aot) + ^alt 2 y (t)f"(a t) + ty (t)f'(aot)a 2 



dt 
(3.25) 



and 



a 2 



1 r 



y2{x)f{a x) + a 1 xy 1 {x)f{aQx) + -a x x y Q {x)f"(a x) 



xy (x) f (a x)dx 



dx 



(3.26) 



Therefore, a second order approximation to the shape of the limit cycle is yo(x) + 
yi(x) + y 2 (x) and a second order approximation to the amplitude of the limit cycle 
is a = a + a\ + a 2 . 

Moreover, at every order in p the appropriate function h(e) is indicated by the 
experiment. 



3.2 Results for the van der Pol equation 



For the van der Pol system, we have f{x) = x 2 — 1, F(x) = x 3 /3 — x, Xq = — 1, 
a = 2 and 



y(x) 







if - 1 < x < -1/2, 



2 \F(2)-F(2x) if -1/2<x<1. 



(3.27) 



The expressions for y n (x) given by the recursive formula (3.20) can be calculated 



with Mathematica. Although they are elementary functions, their length is excessive 



to be reproduced here. From (3.19) we get 



yo(t)f(a t)dt 
-l 64 



96 + 87T, 2 



a 2 Q - 4) = => a = 2. 



(3.28) 



From (3.24) we obtain 



ai(h) 



(3.29) 



35(87r + 9e)' 

Therefore, a first approximation to the amplitude of the limit cycle is 

a(h)~2 r. 3.30 

V ; 35(87r + 9e) V ; 

At it was advanced at the beginning of this section, note that this expression is only 
valid for positive e. Choosing b = —1.3 and c = 4 in (3.13), 

1.3e 



h(e) 



4 + e 2 



(3.31; 



the formula am(e) given in (2.4) is obtained. 



Fig. 3. Approximated shape of the van der Pol limit cycle for e = 5 and h{e) given in 
( |3.34 ). Black line: the experimental limit cycle. Brown, blue and red lines: the curves 



Uo(x), y + y\{x) and y Q + yi(x) + y 2 (x), respectively. 



From (3.26) we get 



a 2 (h) = 

3e/i[19197V3e 2 /z - 5040(8^ - 21h)n + 2e(-22680 v / 3 + /i(54837 + 7700^3vr))] 

19600(9e + 8vr) 2 ' 

(3.32) 

Therefore, a second approximation to the amplitude of the limit cycle is 

a(h) ~ 2 - + 

V ; 35(8vr + 9e) 

3e/i[19197 v / 3e 2 /z - 5040(8^3 - 2lh)n + 2e(-22680 v / 3 + /i(54837 + 7700\/37r))] 

+ 19600(9e + 8tt) 2 ' 

(3.33) 



Choosing b = —0.561 and c = 4 in (3.13) 



, . . 0.561e 

He) = -j^, (3.34) 



the formula a^it) shown in (2.5) is finally obtained. 



As an example, and to conclude this section, we plot in Fig. 3 the form of the van 
der Pol limit cycle when y~o(x), yo{x) + yi(x) and yo(x) + yi(x) + y 2 {x) are used as 
approximated limit cycles. It is not banal to recall here that the reconstruction of 
the shape of limit cycles is not a problem less difficult than that of the calculation 
of their amplitudes. 



4 Conclusions 



In this work, the conjecture (2.8) has been posed. This is a consequence of the 
different formulas here presented for approximating the amplitude a of the van der 
Pol limit cycle. In addition to the well known constant approximation a = 2 that 
generates an error less than 1%, we establish a family of recursive formulas that are 
valid for the whole range of the parameter e. Two of them, ORi(e) and a^i^), have 
been explicitly given. The first one, am, produces an error less than 0.1% and the 
second one, a# 2 , reduces the error to less than 0.05%. Moreover, a R2 is conjectured 
to be an upper bound of a. 

As far as we know, this is the first time where an analytical approximation of the 
amplitude of the van der Pol limit cycle, with validity from the weakly up to the 
strongly nonlinear regime, is proposed. 
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